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Abstract 

The collision of convex bodies is considered for small impact velocity, when plastic deforma- 
tion and fragmentation may be disregarded. In this regime the contact is governed by forces 
according to viscoelastic deformation and by adhesion. The viscoelastic interaction is de- 
scribed by a modified Hertz law, while for the adhesive interactions, the model by Johnson, 
Kendall and Roberts (JKR) is adopted. We solve the general contact problem of convex vis- 
coelastic bodies in quasi-static approximation, which implies that the impact velocity is much 
smaller than the speed of sound in the material and that the viscosity relaxation time is much 
smaller than the duration of a collision. We estimate the threshold impact velocity which dis- 
criminates restitutive and sticking collisions. If the impact velocity is not large as compared 
with the threshold velocity, adhesive interaction becomes important, thus limiting the validity 
of the pure viscoelastic collision model. 



8.1 Introduction 

The large set of phenomena observed in granular systems, ranging from sand and powders on 
Earth to granular gases in planetary rings and protoplanetary discs, is caused by the specific 
particle interaction. Besides elastic forces, common for molecular or atomic materials (solids, 
liquids and gases), colliding granular particles also exert dissipative forces. These forces 
correspond to the dissipation of mechanical energy in the bulk of the grain material as well 
as on their surfaces. The dissipated energy transforms into energy of the internal degrees of 
freedom of the grains, that is, the particles are heated. In many applications, however, the 
increase in temperature of the particle material may be neglected (see, e.g. |6| ). 

The dynamical properties of granular materials depend sensitively on the details of the dis- 
sipative forces acting between contacting grains. Therefore, choosing the appropriate model 
of the dissipative interaction is crucial for an adequate description of these systems. In real 
granular systems the particles may have a complicated non-spherical shape, they may be non- 
uniform and even composed of smaller grains, kept together by adhesion. The particles may 
differ in size, mass and in their material properties. In what follows we consider the contact 
of granular particles under simplifying conditions. We assume that the particles are smooth, 
convex and of uniform material. The latter assumption allows us to describe the particle de- 
formation by continuum mechanics, disregarding their molecular structure. 

It is assumed that particles exert forces on each other exclusively via pairwise mechanical 
contact, i.e., electromagnetic interaction and gravitational attraction are not considered. 
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8.2 Forces Between Granular Particles 
8.2.1 Elastic Forces 

When particles deform each other due to a static (or quasi-static) contact they experience an 
elastic interaction force. Elastic deformation implies that, after separation of the contacting 
particles, they recover their initial shape, i.e., there is no plastic deformation. The stress tensor 
Cgi(r*) describes the i-component of the force, acting on a unit surface which is normal to 
the direction j — {x,y, z}). In the elastic regime the stress is related to the material 
deformation 

where u(r ) is the displacement field at the point r in the deformed body, via the linear relation 



= E i \ Uij(f) - ^S ijU u(f )J + E 2 S llUll (r ) . (8.2) 

Repeated indices are implicitly summed over (Einstein convention). The coefficients E\ and 
E 2 read 

Y Y 

El = 77" — V' E<2 = Tn ' ( 8 - 3 ) 

(1 + v) 3(1 — 2u) 

where Y is the Young modulus and v is the Poisson ratio. Let the pressure P(x, y) act on 
the surface of an elastic semispace, z > 0, leading to a displacement field in the bulk of the 
semispace 1 18 1: 

G lk {x -x',y- y', z) P k (x' , y') dx'dy 1 , (8.4) 

where Git (x, y, z) is the corresponding Green function. For the contact problem addressed 
here we need only the z-component of the displacement on the surface z = 0, that is, we need 
only the component 

°-<-'- »-^^T7-^ 

of the Green function 1 1 8 1 . 

Consider a contact of two convex smooth bodies labeled as 1 and 2. We assume that only 
normal forces, with respect to the contact area, act between the particles. In the contact region 
their surfaces are flat. For the coordinate system centered in the middle of the contact region, 
where x = y = z = 0, the following relation holds true: 



Bxx 2 + B 2 y 2 + u z \(x,y) + u z2 (x,y) = £ , 



(8.6) 
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where u z \ and u z % are respectively the z-components of the displacement in the material of 
the first and of the second bodies on the plane z = 0. The sum of the compressions of both 
bodies in the center of the contact area defines £. The constants B\ and B 2 are related to the 
radii of curvature of the surfaces in contact 1 1 8 1 : 



2 (Bi + B?) =— + — + \ + \ 

ill 1X2 rti 

^ 2 ^ (i?i i?2^ (^-R'i -R2 ) ^O^i -R2 ) O^i R'2 



(8.7) 



Here i?2 and R'±, R' 2 are respectively the principal radii of curvature of the first and the 
second body at the point of contact and ip is the angle between the planes corresponding to 
the curvature radii i?i and R[ . Equations (18.611 . d8.7l > describe the general case of the contact 
between two smooth bodies (see 1 18 1 for details). The physical meaning of (18.61 is easy to 
see for the case of a contact of a soft sphere of a radius R (R\ = R2 = R) with a hard, 
undeformed plane (R[ = R' 2 = 00). In this case B\ = B2 = 1/R, the compressions of the 
sphere and of the plane are respectively u z i(0, 0) = £ and u Z 2 — 0, and the surface of the 
sphere before the deformation is given by z(x,y) = (x 2 + y 2 )/R. Then ( 18 .61 reads in the 
flattened area u z i(x, y) = £ — z(x, y), that is, it gives the condition for a point z(x, y) on the 
body's surface to touch the plane z = 0. 

The displacements u z \ and u Z 2 may be expressed in terms of the normal pressure P z (x, y) 
which acts between the compressed bodies in the plane z = 0. Using (18.4-t and J8.5I > we rewrite 
(IO as 

1 (1 -1 , 1 - 4\ rr p *(x'>y') dx > dy > = { _ BlX * B ^ , (8 . 8) 



7T V Fl Y 2 



where r = y/(x — x') 2 + (y — y') 2 and integration is performed over the contact area. Equa- 
tion (18.81 is an integral equation for the unknown function P z (x,y). We compare this equation 
with the mathematical identity 1 1 8 1 



dx'dy' / x' 2 y' 2 nab f \ x 2 y 2 



b 2 2 



a 2 +t b 2 +t 



dt 

(8.9) 



^{a 2 +t){b 2 +t)t 



where integration is performed over the elliptical area x' 2 /a 2 + y' 2 /b 2 — 1. The left-hand 
sides of both equations contain integrals of the same type, while the right-hand sides contain 
quadratic forms of the same type. Therefore, the contact area is an ellipse with the semi-axes 
a and b and the pressure is of the form P z (x, y) = const-\/l — x 2 ja? — y 2 /b 2 . The constant 
may be found from the total elastic force F c i acting between the bodies. Integrating P z (x, y) 
over the contact area we obtain 



We substitute (18. 10l > into J8.8I > and replace the double integration over the contact area by 
integration over the variable t, according to the identity (18 .9b . Thus, we obtain an equation 
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containing terms proportional to x 2 , y 2 and a constant. Equating the corresponding coeffi- 
cients we obtain 



Bi = 
B 2 = 



F cl D 



dt 



F el D N(x) 



TT 

F e lD 
TT 

FelD 
TT 



o ^{a 2 +t){b 2 + t)t it 
dt 



FeiD M(x) 



o (a 2 + t)^/(a 2 +t)(b 2 + t)t 7T a 2 b 

dt F el DM(l/x) 



where 



D = 



o {b 2 +t)^{a 2 + t)(b 2 +t)t * 



1-^2 

Y 2 



ib 2 



(8.11) 
(8.12) 
(8.13) 

(8.14) 



and x = a 2 /b 2 is the ratio of the contact ellipse semi-axes. In ( 18.1 H -l l8T3l we introduce the 
short-hand notations 1 



N(x) 
Mix) 



dt 



+Xt)(l +t)t 

dt 



(8.15) 
(8.16) 



/O (l+t)y/(l+t)(l+Xt)t 

From these relations will follow the size of the contact area, a, b, and the compression £ as 
functions of the elastic force F c \ and the geometrical coefficients B\ and B 2 - 

The dependence of the force F c \ on the compression £ may be obtained from scaling 
arguments. If we rescale a 2 — ► aa 2 , b 2 — > ab 2 , £ — * a£ and F c i — > a 3 ' 2 -F c i, with a constant, 
Eqs. J8 . 1 1 i — <l8 . 1 3i remain unchanged. That is, when £ changes by the factor a, the semi-axis 
a and b change by the factor a 1 / 2 and the force by the factor a 3 / 2 , i.e., a <~ t; 1 / 2 , b ~ t; 1 / 2 
and 



Fei = const £ 



3/2 



(8.17) 



The dependence 18.171 holds true for all smooth convex bodies in contact. To find the constant 
in 18. 171 we divide 18. 131 by 18 . 1 21 and obtain the transcendental equation 



Eh = y/xM (1/x) 
Si ~ M(x) 



(8.18) 



for the ratio of semi-axes x. Let xo be the root of Eq. (18. 1 8i . then a 2 = xob 2 and we obtain 

F c iD N{x ) 



B 1 = 



TT b 
FelD M(x ) 
TT Xgb 3 



(8.19) 
(8.20) 



'The function N(x) and M(x) may be expressed as a combination of the Jacobian elliptic functions E(x) and 
K(x) 
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where N(xq) and M(xo) are pure numbers. Equations ( I8.19> . J8.20i allow us to find the 
semi-axes b and the elastic force F c \ as functions of the compression £. Hence we obtain the 
force, i.e., we get the constant in ( I8.17t : 

1/2 



7T / M(X ) \ 3/2 = 3/ 

D \Bix N(xa)J ? ° ? 



*d = ( TT~~"ir77" — \ ) t"* = Co? 12 • (8.21) 



For the special case of contacting spheres (a = 6), the constants B\ and B2 read 

Si=B H(i + £Hi^ (8 - 22) 

In this case x = 1, N(l) = n, and M (1) = ?t/2, leading to the solution of J8.19L (18.20b : 

a 2 = R cS £ (8.23) 

f cI = P? 3/2 ; p^^-^VR^, (8.24) 

where we use the definition (18.141 of the constant D. This contact problem was solved by 
Heinrich Hertz in 1882 [ 14 1. It describes the force between elastic particles. For inelastically 
deforming particles it describes the repulsive force in the static case. 



8.2.2 Viscous Forces 

When the contacting particles move with respect to each other, i.e., the deformation changes 
with time, an additional dissipative force arises, which acts in the opposite direction to the 
relative particle motion. The dissipative processes occurring in the bulk of the body cause 
a viscous contribution to the stress tensor. For small deformation the respective component 
of the stress tensor is proportional to the deformation rate ii.y(r), according to the general 
relation 181: 



<r % j s J?,t)=E l ]drnl> 1 {t-T) 




+ E 2j dTijj 2 (t-T)8ijU U (f, r) , 


(8.25) 



where the (dimensionless) functions ipi(t) and ifaft) are the relaxation functions for the dis- 
tortion deformation and ^2 (i) for the dilatation deformation. 

In many important applications the viscous stress tensor may be simplified significantly. 
If the relative velocity of the colliding bodies is much smaller than the speed of sound in the 
particle material and if the characteristic relaxation times of the dissipative processes r v ; s xj2 
are much smaller than the duration of the collision t r , 



,1/2 



(8.26) 



the viscous constants 771 and 772 may be used instead of the functions ipi(t) and ifaif). Thus 

(8.27) 



Vl/2 = E 1/2 T visA /2 = E±/2 / lp 1 / 2 (T)dT 

Jo 
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and the dissipative stress tensor reads (see 1 8 1 for details) 

1 



Uij(r,r) - -5ijiiu(r,T) 



(8.28) 



It may be also shown that the above conditions are equivalent to the assumption of quasi- 
static deformation |8 7|. When the material is deformed quasi-statically, the displacement 
field u(r ) in the particles coincides with that for the static case u c \(f ), which is the solution 
of the elastic contact problem. The field u c \(r ), in its turn, is completely determined by the 
compression £, which varies with time during the collision, i.e., u c \ — u c \(r, £). Therefore, 
the corresponding displacement rate may be approximated as 



• d 

u{r,t) ~ £-^-Uei(r,£) 
and the dissipative stress tensor reads, respectively 



l 



(8.29) 



(8.30) 



From (18.301 and J 8 . 2 1 follows the relation between the elastic and dissipative stress tensors 
within the quasi-static approximation, 



(8.31) 



where we emphasize that the expression for the dissipative tensor may be obtained from the 
corresponding expression for the elastic tensor after substituting the elastic constants by the 
relative viscous constants, and application of the operator £d/d£. 

The component er^f of the elastic stress is equal to the normal pressure P z at the plane 
z = of the elastic problem, Eq. (18. lOt 



* z J(x,y,0) =25!— * 
oz 



du x 
dx 



du y 
dy 



duz 

dz 



. 3F cl 
2 nab 



(8.32) 



Now we compute the total dissipative force acting between the bodies. Instead of a direct 
computation of the dissipative stress tensor, we employ the method proposed in |E]EQ : We 
transform the coordinate axes as 



ax 



y = ay 



(8.33) 



with 







(E 2 + 




\m + 




\E 2 - 


kEj 



(3 = 



(V2 



hi) 



a(E 2 - ii5i) 



b = ab' 



(8.34) 
(8.35) 
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and perform the transformations 



du z 
dz 

= P 
= P 



Ei 



m - 

du- 
~dz~' 



m 

3 



du x du y 
dx dy 
Ex 
3 



duz 
dz 



( du x 
{'dx' 



3F< 



el 



2 ira'b' 



,,/2 



„/2 



1 



771 - 772 = P<* 



8Uy 

H 2 

3F e i 
2 7ra6 



9u 2 
~dz' 



(8.36) 



^ a; 2 ?/ 2 



b 2 



Applying the operator £,d/d£ to the last expression on the right-hand side we obtain the dis- 
sipative stress tensor. Subsequent integration over the contact area yields, finally, the total 
dissipative force acting between the bodies: 



where 



2 R 1 (3i7a 
a p = — 



3 (3?7 2 + 2 m 



(1 



(l-2i/) 



y^ 2 



(8.37) 



(8.38) 



Using the scaling relations for the elastic force, Eq. J8.17> . and for the semi-axes of the contact 
ellipse, we obtain 



dF cl 



3£k 
2 £ 



da 



la 

H 



db 



lb 
2£ 



(8.39) 



Then from J8.36i and (18.2 II . the distribution of the dissipative pressure in the contact area 
may be found: 



Air ab 



-1/2 



(8.40) 



where the constant Co is defined in J8.2H . 

We wish to stress that, to derive the above expressions, we assumed only that the surfaces 
of the two bodies in the vicinity of the contact point before the deformation, are described 
by the quadratic forms Z\ — nf^XiXj and z 2 = K^XiXj — x,y, z), where k^ 2 '' are 
symmetric tensors 1 18 1. Therefore, the relations obtained are valid for a contact of arbitrarily 
shaped convex bodies. For spherical particles of identical material, J8.37t and J8.24l > yield 



0171 



(8.41) 



with p as defined in J8.24I . Hence, the total force acting between viscoelastic spheres takes 
the simple form |8][7J 



F = p[t 



r3/2 



(8.42) 
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The range of validity of J8.42> for the viscoelastic force is determined by the quasi-static 
approximation. The impact velocity must be significantly smaller than the speed of sound. On 
the other hand, the impact velocity must not be too small in order to neglect adhesion. We 
also neglect plastic deformation in the material. 

8.2.3 Adhesion of Contacting Particles 
Models of Adhesive Interaction 

The Hertz theory has been derived for the contact of non-adhesive particles. Adhesion be- 
comes important when the distance of the particle surfaces approaches the range of molecular 
forces. Johnson, Kendall and Roberts (JKR) 11171 extended the Hertz theory by taking into 
account adhesion in the fiat contact region. They show that the contact area is enlarged by 
the action of the adhesive force. Therefore, they introduced an apparent Hertz load Fjj which 
would cause this enlarged contact area. To simplify the notation, we consider the contact 
of identical spheres. The contact area is then a circle of radius a, which corresponds to the 
compression for the Hertz load Fh- In reality, however, this contact radius occurs at the 
compression £ which is smaller than In the JKR theory it is assumed that the difference 
between the Hertz compression and the actual one, £, may be attributed to the additional 
stress 



which is the solution of the classical Boussinesq problem 1 32 1 . This distribution of the normal 
surface traction gives rise to a constant displacement over a circular region of an elastic body. 
The displacement £g corresponding to the contact radius a and the total load Fb are related by 



where the constant D is defined in (18.14V 

The value of Fb < mimics the additional surface forces, such that the pressure is positive 
(compressive) in the center of the contact area, while it is negative (tensile) near the boundary 
1171 . Hence, the shape of the body is determined by the action of two effective forces Fh and 
Fb- The total force between the particles is their difference, F = Fh — Fb- Johnson et al. 
assumed that the elastic energy stored in the deformed spheres may be found as a difference 
of the elastic energy corresponding to the Hertz force Fjj and that due to the force Fb 1171 . 
Using 




(8.43) 




(8.44) 



U s = —irja 2 



(8.45) 



for the surface energy, where 7 > is twice the surface free energy per unit area of the solid 
in vacuum or gas, and minimizing the total energy, we find 11171 
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a 3 = \DR |f + '^ttjR+ \I3itjRF+ [ ^ttjR ) | (8.47) 
and also the compression 



The first term in (I8.48i is the Hertz compression which coincides with J8.23l > for R cS — > 
R/2. Equation (18.471 may be solved to express the total force as a function of the contact 
radius: 



2a A /6tt7 3/2 



For vanishing applied load the contact radius ao is finite: 

<4 = ^DitjR 2 . (8.50) 

For negative applied load the contact radius decreases and the condition for a real solution of 
J8.47t yields the maximal negative force which the adhesion forces can resist, 



F sep = -^7T7./,\ (8.51) 



corresponding to the contact radius 



For a larger (in the absolute value) negative force, the spheres separate. For spheres of dissim- 
ilar radii, in (l8~47l - d8~52l R should be substituted by 2R eS . 

Another approach to the problem of the adhesive contact was developed by Derjaguin, 
Muller and Toporov (DMT). They assumed that the Hertz profile of the pressure distribution 
on the surface stays unaffected by adhesion and obtained the pull-off force F sop = —2ttjR bS 
(9). The assumption of the Hertz profile allows one to avoid the singularities of the pressure 
distribution (I8.43> on the boundary of the contact zone. Since the experimental measurement 
of 7 is problematic, it is not possible to check the validity of the JKR and DMT theories, i.e., 
to resolve their disagreement. 

In later studies 1201 1211 a more accurate theoretical analysis has been performed. The 
elastic equations have been solved numerically for a simplified microscopic model of adhe- 
sive surfaces with Lennard-Jones interaction. Within this microscopic approach, the relative 
accuracy of different theories has been estimated for a wide range of model parameters. It 
was found that the DMT theory is valid for small adhesion and for small, hard particles. JKR 
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theory is more reliable for large, soft particles with large adhesion forces, which, however, 
should be short-ranged. 

In 1 2 1 the Lennard-Jones continuum model of solids was studied. The adhesive forces 
between the surfaces then read 



— - 1 



(8.53) 



Here P s (h) describes the forces acting per unit area between the surfaces, h = h(r) is the 
actual microscopic distance between them. H is the Hamaker constant, characterizing the van 
der Waals attraction of the particles in a gas or vacuum and zo is the equilibrium separations 
of the surfaces. The surface energy in this model is defined by 

7 = t^-t • (8-54) 

It was observed in 1 2 1 that the accuracies of different theories vary depending of the value of 
the Tabor parameter /i, 1 29 1 



,i 3 / 2 = - 

M 3 



jDyJlFB/z*. (8.55) 

In agreement with |20 21 1 it has been shown |2 1 that small values of /i (small hard particles 
with low surface energies) favor the DMT theory (/i < 10~ 2 ) while for \i ~ 1-10 the JKR 
theory proves to be rather more accurate. Both JKR and DMT fail for large \i when the 
strong adhesion is combined with the soft material of the contacting bodies. In this limit, the 
surfaces jump into contact, which corresponds to a spontaneous non-equilibrium transition 
(see e.g. |27|). Similar analysis has been performed later [IX], where the author concluded 
that the DTM theory generally fails, both in original and corrected forms. One of the main 
conclusions of El ll II is that the JKR theory, albeit simple, gives relatively accurate predictions 
for basic quantities in the range of its validity (/i ~ 1-10). 

Among the theories developed to cover the DMT-JKR transition regimes |20 2Tl ll lll29l 
1161 1191 the theory by Maugis |19| is the most frequently used. It is based on a simplified 
model of adhesive forces. The adhesive force of constant intensity Pjj is extended over a 
fixed distance fi£> above the surface, yielding the surface tension 7 = Pjj hp . The description 
of a contact in this model is based on two coupled analytical equations which are to be solved 
numerically. The recently developed double-Hertz model II 2111 31 constructs the solution for 
the adhesive contact as a sum of two Hertzian solutions, which make the theory analytically 
more tractable than the Maugis model. Combining, in the adopted manner, the successful 
assumptions of the JKR and the modified DMT model, a generalized analytical theory for the 
adhesive contact has been proposed 1 26 1. 

In what follows we assume that the parameters of our system belong to the range of validity 
of the JKR model, fj, ~ 1 — 10, and will use this simple analytical theory to describe the 
adhesive contacts between spheres. Moreover, we assume that the adhesive force is small. 
To estimate the influence of the adhesive force, we approximate £ as 2a 2 /R in J8.48> and 
substitute it into J8.49I I to obtain (see also 1281 ~). 

F « P e /2 - V&H/D (R cS ) 3/4 e /A ■ (8.56) 
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Viscoelasticity in Adhesive Interactions 

The adhesive forces between particles cause the additional deformation in the contacting bod- 
ies as compared to a pure Hertzian deformation, hence in the corresponding dynamical prob- 
lem an additional deformation rate arises. Therefore, the dissipative forces must have an 
additional component attributed to the adhesive interactions. The adhesive contact of vis- 
coelastic spheres has been studied numerically in 1281 1131 . In the quasi-static condition 
for the colliding viscoelastic adhesive spheres was used and an analytical expression for the 
interaction force has been derived for the JKR model. Similar to the case for non-adhesive 
particles, it was assumed that in the quasi-static approximation, the deformation field may 
be parameterized by the value of the compression £. (Note that this assumption neglects the 
possible hysteresis which can happen for the negative total force Q). Performing the same 
transformation which lead to the expression J8.37t for the case of non-adhesive contact, and 
using the approximation J8.56i we obtain the estimate for the dissipative forces |5 1 

F dis = \ApiJl+ iBy/^/D (i? cff ) 3/4 Cr 1/4 (8.57) 

B ee aP = faz^jL . (8.58) 
h 3(l + z/)(l-2i/) 

Note the singularity in the second term of J8.57l > at £ = 1 . It is attributed to the quasi-static 
approximation for JKR theory and physically reflects the fact that the adhesive particles can 
jump into contact 1127 1 with the discontinuous change of the compression £. Consider now 
how the above forces determine the particle dynamics. 



8.3 Collision of Granular Particles 

8.3.1 Coefficient of Restitution 

Based on the particle interaction forces discussed so far, we turn now to the description of 
the particle collisions. It is assumed that the colliding particles do not exchange tangential 
forces 3 , hence, only normal motion is considered. Let the particles be spheres of the same 
material, which start to collide at time t = at relative normal velocity g (impact rate). The 
time-dependent compression then reads 

£(t) =Ri + Rj - \n(t) - fj{t)\ , (8.59) 

where rj(t) and fj(t) are positions of the particle centers at time t (see Figure lOt . The 
relative normal motion of particles at a collision is equivalent to the motion of a point particle 
with the effective mass 



2 This is a weak or integrable singularity, that is Jq£~ 1//4 c££ ~ e 3 / 4 — » for e — » 0. Hence for practical 
application of 18.571 one can use § > e, where e may be very small but a finite number. 
3 See 1 6 for a discussion of the consistency of this assumption. 
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Figure 8.1: Head-on collision of identical spheres. The time-dependent state is characterized by the 
compression £(£) = 2R — \ fi(t) — r2{t)\ and the compression rate f(t) = Vi(t) — V2(t). 

For the moment let us neglect adhesion and consider the collision of viscoelastic particles 
interacting via the force (I8.42> . The equation of motion and the initial conditions read 

e+^i f^ + fVee) -o, m= g , e(o) = o. (8.6D 

When granular particles collide, part of the energy of the relative motion is dissipated. The 
coefficient of (normal) restitution quantifies this phenomenon: 

e = -i(t e )/t(0) = -i(t e )/g, (8.62) 

where £ (0) = g is the pre-collision relative velocity and t c is the duration of the collision. 
In general, e is a function of the impact velocity. It can be obtained by integrating J8.6H 
numerically OBOE) or analytically l24l . 

8.3.2 Dimensional Analysis 

The analytical solution 1 24 1 requires considerable efforts, here we give a simplified derivation 
which is based on a dimensional analysis of the equation of motion J8.6H [23 1. This method 
was employed before 1311 to prove that the frequent assumption e =const. is inconsistent 
with mechanics of materials. For the dimensional analysis, the elastic and dissipative forces 
are represented in the more general form 

F cl = m eS Dt C , F diss = m cfi D 2 , (8.63) 

with D1/2 being material parameters. With these notations the equation of motion for collid- 
ing particles reads 

^ + D 1 C + D 2 Ci P = 0, £(0)=g, C(0)=0. (8.64) 

For the case of pure elastic deformation (D 2 = 0) the maximal compression £o is obtained by 
equating the initial kinetic energy, m eS g 2 /2 and the elastic energy m cS Di^ +1 / (a + 1): 

/n- + 1\ l/(l+Q) 
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We chose £o as the characteristic length of the problem. The time needed to cover the dis- 
tance £o when traveling at velocity g defines the characteristic time: tq = £o/.9- Thus, the 
dimensionless variables read 



In dimensionless form, J8.6H reads 



i £0 



l + a 



r = o. 



e=(Wfl 2 )e- 



1(0) = 1 , 



1(0) = o 



(8.66) 



(8.67) 



with 



x = x(g) = D 2 



1 + cA 
2D X ) 



(l+ 7 )/(l+a) 



9 



2( 7 -a)/(l+a)+/3 



(8.68) 



None of the terms in J8.67t depends either on material properties or on impact velocity, except 
for x. Therefore, if the motion of the particles depends on material properties and on impact 
velocity, it may depend only via n, i.e., in the combination of the parameters as given by 
J8.68i . Hence, any function of the impact velocity, including the coefficient of restitution 
must be of the form e(g) = e[x(g)]. A similar result for e — > 0, (3 = 1 and a — 3/2 has been 
obtained in 1 10 1. 

Hence, if the coefficient of restitution does not depend on the impact velocity g, it is 
implied that 



2( 7 -a)+/3 (l + a) = 0. 



(8.69) 



For small £ a linear dependence of the dissipative force on the velocity seems to be realistic, 
i.e., (3 = 1. Then e =const. holds true for the following cases: 



• For the linear elastic force F c \ oc £, (i.e. a 
dashpot force Fdi s oc £, (7 = 0). 



1) condition J8.69i implies the linear 



• For the Hertz law for 3D-spheres J8.24> . (i.e. a — 3/2), condition (18.691 requires F^is oc 
CC 1 ^ 4 , (7 = |f). As far as we can see there is no physical argument to justify this 
functional form of the dissipative force. 

Therefore, we conclude that the assumption e =const. is in agreement with mechanics of 
materials only in the case of (quasi-)one-dimensional systems. For three-dimensional spheres 
it disagrees with basic mechanical laws. 

For viscoelastic spheres, according to J8.42t . the coefficients are a = 3/2, /3 = 1, and 
7 = 1/2. From J8.68l > it follows that 



and, therefore, 



.4 



m J 



TO / 



(8.70) 



(8.71) 
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If we assume that the function e{g) is sufficiently smooth and can be expanded into a Taylor 
series, and with e(0) = 1, for small impact velocity the coefficient of restitution reads 

e = l- dAnWg 1 '* + C 2 A 2 K 4 / 5 g 2 ^ T ■ ■ ■ (8.72) 

where 



(8.73) 



m eff (1 - v 2 ) ' 

The coefficients C\, Ca, . . . are pure numbers which are given analytically in 1241 . Here we 
give a simple derivation of these coefficients (which is correct for C\ and C2 and approxi- 
mately correct for C3 and C4, using the method proposed in 1 23 1). 



8.3.3 Coefficient of Restitution for Spheres 
Small Inelasticity Expansion 

Using d/d£ = £d/d£ the equation of motion for a collision adopts the form 

K^ + 5^ /9 ) = "*^ = ^. * (0) = ' ^ (0) = 1 ' (8 ' 74) 



where we introduce the mechanical energy 



(8.75) 



The first stage of the collision starts at £ = and ends in the turning point of maximal 
compression £0. During the second stage, the particles return to £ = 0. The energy dissipation 
during the first stage is given by 

so Jf rio ■ r- 

-jdd = -K t^dt. (8.76) 
d£ Jo 

For the evaluation of the right-hand side of J8.76> . the dependence £ = £(£) is needed. In 
the case of an elastic collision where the maximal compression is £0 = 1 (according to the 
definition of the dimensionless variables) from energy conservation, it follows that 



£(£) = yi-P' 2 , < 8 - 77 ) 

i.e., £ vanishes at the turning point £ = 1. For inelastic collisions £0 < 1, therefore, 



Using (I8.78i the integration in (I8.76> may be performed yielding 

\iv' 2 -\ = -^^ 12 (8.79) 
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where we take into account 

E(io) = \H' 2 , E(0) = || 2 (0) = \ (8.80) 

and introduce the constant 

Let us define the inverse collision, the collision that starts with velocity e g and ends with 
velocity g. During the inverse collision, the system gains energy. The maximal compression £o 
is naturally the same for both collisions, since the inverse collision equals the direct collision, 
except for the fact that time runs in the reverse direction, hence, 

^jf-=+4\ft, 1(0) =e, |(0) = 0. (8.82) 

This suggests an approximative relation for the inverse collision, 



etewv 1 -^") 872 ' (8 - 83) 

with the additional pre-factor e, which is the initial velocity for the inverse collision. 

Integration of the energy gain for the first phase of the inverse collision (which equals, up 
to its sign, the energy loss in the second phase of the direct collision |24|) may be performed 
just in the same way as for the direct collision, yielding 

5& 5/2 - ^=+SHbil / \ (8.84) 



it with (18.841 . the maximal compression is e — £ ■ Substituting this into J8.79i we arrive at 



where again E(i ) = |q /2 /2 and E{0) = e 2 /2 is used. Multiplying ( 18791 by e and summing 
it with (18.841 . the maximal compression is e 
an equation for the coefficient of restitution 

£ + 2xbe 3/5 = 1. (8.85) 

The formal solution to this equation may be written as a continuous fraction (which does not 
diverge in the limit g — > oo): 

e" 1 = 1 + 2x6(1 + 2x6(1 + • • • ) 2/5 • • • ) 2/5 (8.86) 

For practical applications the series expansion of e in terms of h is more appropriate. We 
return to dimensional units and define the characteristic velocity g* such that 



1 / 9 x 



26 \g 



(8.87) 
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with b being defined in J8.8U . Using, moreover, the definition (I8.68> together with J8.42I . 
which provides the values of D\ and D2, the characteristic velocity reads 

, * r i/ 5 _ ^ r(3/5) (3 \ 

(3) - 21/552/5 r (21/10) U J ■ 

With this new notation the coefficient of restitution adopts a simple form: 

/ \ 1/5 / \ 2/5 / x 3/5 / x 4/5 

with ai = 1, a 2 = 3/5, a 3 = 6/25 = 0.24, a 4 = 7/125 = 0.056. Rigorous but elaborated 
calculations |24| show that, while the coefficients a\ and are exact, the correct coefficients 
03 and 04 are: 03 ~ 0.315 and 04 » 0.161. The coefficients Cj of the expansion (18.721 can 
be obtained via 

Q = ai C{ = ai (g*y l/5 . (8.90) 
In particular, 

y¥ r(3/5) _ 3 2 

6 i- 21/552/5 r (21/10)' 62 5 1 (8 - 91) 

and respectively, C3 » —0.483, C4 » 0.285. The convergence of the series is rather slow, 
and accurate results can be expected only for small enough g/g*. 

Let us briefly mention a complication of the quasi-static approximation (QSA). During the 
expansion phase it may happen that the repulsive force according to (18.421 becomes negative, 
i.e., seemingly the particles attract each other. For the interaction of non-cohesive particles we 
had, however, excluded attractive forces. This is an artefact, since in reality the particles lose 
contact already, before completely recovering their spherical shape, i.e., before £ = (see 
1 22 1 for a detailed explanation of this problem). This effect, however, is not in agreement with 
the QSA. Obviously, 18.421 which is a result of the QSA, derived in [8|, is not appropriate 
to describe the very end of the particle contact. Taking this effect into account we obtain a 
larger coefficient of restitution as compared with the presented computation 1251 . For small 
dissipation, the correction is rather small. This small correction is neglected here. 

Pade Approximation 

For practical applications, such as molecular dynamics simulations, the expansion ( I8.89i is of 
limited value, since it diverges for large impact velocities, g — ► 00. It is possible, however, 
to construct a Pade approximant for e, based on the above coefficients, which reveals the 
correct limits, e(0) = 1 and e(oo) = 0. The dependence of e(g) is expected to be a smooth 
monotonically decreasing function, which suggests that the order of the numerator must be 
smaller than the order of the denominator. The 1-4 Pade-approximant 

£ = ! + <*! (g/fl*) V5 (8 92) 

1 + d 2 (g/g*) 1/5 + d 3 (g/g*f 5 + d 4 (g/g*f 5 + d 5 (g/g*f 5 
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Figure 8.2: Dependence of the coefficient of normal restitution on the impact velocity for ice particles. 
The dashed line is experimental |4|, the solid line is the Pade-approximation J8.92I with the constants 
given by 18.931 and with the characteristic velocity for ice g* = 0.32 cm s -1 . 

satisfies these conditions. Standard analysis (e.g. [3 1) yields the coefficients dh in terms of the 



coefficients eife 

d Q = a 4 - 2a 3 - aj + 3a 2 - 1 (8.93) 

di = [1 - a 2 + a 3 - 2a 4 + (a 2 - l)(3a 2 - 2a 3 )] /d ~ 2.583 

d 2 = [(«3 - a 2)(l - 2a 2 ) - a 4 ] /d w 3.583 

d 3 = [«3 + a 2 (a 2 - 1) - a4(a 2 + 1)] /d ~ 2.983 

rf 4 = [a 4 (a3 - 1) + (a 3 - a 2 )( a 2 - 2a 3 )] /d ~ 1.148 



rf 5 = [2(03 - a 2 )(a 4 - a 2 a 3 ) - (a 4 - o 2 ) 2 - a 3 (a 3 - a 2 )] /d ~ 0.326 

Using the characteristic velocity g* = 0.32 cm s _1 for ice at very low temperature as a fitting 
parameter, we compare the theoretical prediction of e(g), given by J8.92> . with the exper- 
imental results 1 4 1, see Figure liOl The discrepancy with the experimental data at small g 
follows from the fact that the extrapolation expression, e = 0.32/g 234 used by |4| to fit the 
experimental data has an unphysical divergence at g — > and does not imply the failure of 
the theory for this region. The scattering of the experimental data presented by |4) is large 
for small impact velocity, according to experimental complications, therefore the fit formula 
of 1 4 1 cannot be expected to be accurate enough for velocities that are too small. For very 
high velocities the effects, such as brittle failure, fracture and others, may contribute to the 
dissipation, so that the mechanism of the viscoelastic losses could not be the primary one. 
In the region of very small velocity, other interactions than viscoelastic ones, e.g., adhesive 
interactions, may be important. 

8.3.4 Coefficient of Restitution for Adhesive Collisions 

For very small velocities, when the kinetic energy of the relative motion of colliding particles 
is comparable with the surface interaction energy at the contact, the adhesive forces play an 
important role in collision dynamics - they may change the coefficient of restitution quali- 
tatively. Indeed, as described above, adhesive particles in contact are compressed even for 



Pade approximation 

Experiment by Bridges et al. (1984) 
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vanishing external load, i.e., a tensile force must be applied to separate the particles. There- 
fore, at the second stage of the collision, the separating particles must overcome a barrier due 
to the attractive interaction, which keeps them together. The work against this tensile force 
reduces the kinetic energy of the relative motion after the collision, that is, it reduces the ef- 
fective coefficient of restitution. For small impact velocity the kinetic energy of the relative 
motion may be too small to overcome the attractive barrier, i.e., the particles stick together 
after the collision, corresponding to e = 0. 

From these arguments it follows that the description of particle collisions by pure vis- 
coelastic interaction has a limited range of validity, not only for large impact rate when plastic 
deformation becomes important, but also for small impact rata due to adhesion. A simplified 
analysis of adhesive collisions is presented in [5 1 to estimate the influence of adhesive forces 
on the coefficient of restitution. It allows to estimate the range of validity of the viscoelastic 
collision model. 

We assume that the JKR theory is adequate for the given system parameters. We also 
assume that the adhesion is small and that the adhesive interactions may be neglected when 
the force between the particles is purely repulsive. Hence, we take into account the influence 
of adhesive interaction only when the total force is attractive, that is when the force is mainly 
determined by adhesion. This happens at the very end of the collision. We also neglect the 
additional dissipative forces, which arise due to the adhesive interaction and assume that all 
dissipation during the collision may be attributed to the viscoelastic interactions. 

At the second stage of a collision, when the particles move away from each other they 
pass the point where the contact area is cto and the total force vanishes. As the particles move 
away further, the force becomes negative, until it reaches, at a = a scp , the maximum negative 
value F — F scp , Eq. (I8.5H . At this point the contact of the particles is terminated and they 
separate. According to our assumption, the work of the tensile force which acts against the 
particles, separation reads 



Using J8.491 for the total force F(a), J8.481 for the compression, which allows one to obtain 
d£/da, and J8.50i . J8.52l > for ao and a scp , we obtain the work of the tensile forces, 



Close to the end of the collision, just before the tensile forces start to act, the relative velocity 
is g' = eg. The final velocity g", when the particles completely separate from each other, may 
be found from the conservation of energy: 




(8.94) 




(8.95) 



with the constant 




(8.96) 




m cS (g"f = W ( 



(8.97) 
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From the latter equation we obtain the coefficient of restitution for the adhesive collision, e a( j, 



g" y £ 2 (g). 9 2 - 2Wo/m° ff 

£ad(ff = — = 

9 9 



(8.98) 



where e(g) is the coefficient of restitution without the adhesive interaction. Hence we obtain 
the condition for the validity of the viscoelastic collision model, 



e{g)9 > 



,cff ' 



(8.99) 



The threshold impact velocity g st which separates the restitutive (g > g st ) from the sticking 
(g < g s t) collisions, may be obtained from the solution of the equation 



\m°*e\g)g* = W . 



(8.100) 



Using J8.72i we obtain for viscoelastic spheres, in the leading-order approximation, with re- 
spect to the small dissipative parameter A: 



g s t 



2Wo 



1 + CiAk 2 ' 5 



2W 



1/10 



(8.101) 



For head-on collisions (vanishing tangential component of the impact velocity) the colliding 
particles stick together if g < g at . In this case, after the collision, the particles form a joint 
particle of mass mi + rri2. 



8.4 Conclusion 

We have considered the collision of particles in granular matter with respect to viscoelastic and 
adhesive interaction. Thus, the elastic contribution due to the classical Hertz theory is comple- 
mented by the simplest model for dissipative material deformation, where the viscous stress 
is linearly related to the strain rate. Moreover, quasi-static approximation was assumed, i.e., 
the impact velocity is much smaller than the speed of sound in the material and the viscosity 
relaxation time is much smaller than the duration of the collision. Using these approxima- 
tions, we obtained the general solution for the contact problem for convex viscoelastic bodies. 
The validity of this model is violated for large impact velocity due to plastic deformations 
and also for very small impact velocity due to surface forces. We have discussed the available 
models of adhesive interaction. For the model by Johnson, Kendall and Roberts 1 17 1 which 
has been shown to be accurate in a range of parameters of practical interest, the additional 
dissipation arising due to adhesive forces have been estimated. From the comparison of the 
force contribution due to pure viscoelastic interaction and the contribution due to adhesion, 
we have estimated the range of validity of the viscoelastic model. For head-on collisions we 
have also estimated the marginal value of the impact velocity, which discriminates restitutive 
and sticking collisions. 
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